%Estimate the policy function based on the bootstrapped sample and forward
%simulate future states, including decreasing sunk costs



tN=2009;
t0=2006;
tbar=tN-t0;


%&&& Import the inputs generated by STATA
choice=1:13;
ns=numel(choice);
beta=0.95;


textFileNameY=sprintf('matY.csv');  %the input file
textFileNameX=sprintf('matX.csv');   %the input file
textFileNameZ=sprintf('matZ.csv');   %the input file

tempY=dlmread(textFileNameY);  %matrix: [ahaid year hsa adj_vnd prechs sysid]
tempY=sortrows(tempY,[1,2]);
tempX=dlmread(textFileNameX);  %matrix: [ahaid year hsa bdtot prechs nprofit profit teach hhi totpop65 perc_mcr perc_mcd mktcat2 mktcat3 mktcat4 yrsince sysid]
tempX=sortrows(tempX,[1,2]);
tempZ=dlmread(textFileNameZ);  %matrix: [ahaid year vid hsa vmshbd1 vmshbd2 vmshbd3 vmshbd4 vlagbd1 vlagbd2 vlagbd3 vlagbd4 newvlag newvsysh below50cat below75cat below100cat yrsince_dom prechs sysid]
tempZ=sortrows(tempZ,[1,2,3]);


%&&& Prepare inputs for the four types: SANA, SAAD, AFNA, AFAD
textFileName_mkt=sprintf('smplmkt.csv'); %!!!NOT yet done: Need to re-number for the sample markets used!!!
thisboot=sort(dlmread(textFileName_mkt))';  %row vector
reorder_mkt=1:numel(thisboot); %row vector
baseY=[];
baseX=[];
baseZ=[];
for i=1:numel(thisboot)
    % construct new Y
    onemkt_Y=tempY(tempY(:,3)==thisboot(i),:);
    onemkt_Y(:,3)=reorder_mkt(i);
    baseY=[baseY;onemkt_Y]; 
    % construct new X
    onemkt_X=tempX(tempX(:,3)==thisboot(i),:);
    onemkt_X(:,3)=reorder_mkt(i);
    baseX=[baseX;onemkt_X];    
    % construct new Z
    onemkt_Z=tempZ(tempZ(:,4)==thisboot(i),:);
    onemkt_Z(:,4)=reorder_mkt(i);
    baseZ=[baseZ;onemkt_Z];  
end

% ALSO need to include affiliated hospitals' brothers outside the selected
% markets in the bootstrap sample; also need to reorder mkt ID!!!
sysid=unique(baseY(:,end));
sysid(sysid==0)=[];
baseY_sys=[];
baseX_sys=[];
baseZ_sys=[];
for i=1:numel(sysid)
    % construct Y for affiliated hospitals' brothers outside the selected
    % markets; ALSO need to reorder sysid!!!
    onesys_Y=tempY(tempY(:,end)==sysid(i),:);
    incmkt=onesys_Y(:,3);  %columnn vector
    exclu=max((~(incmkt(:,ones(numel(thisboot),1))-thisboot(ones(size(onesys_Y,1),1),:))),[],2);  %only include bros OUTSIDE the selected markets
    onesys_Y(exclu==1,:)=[];
    %vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
    if isempty(onesys_Y)==0
       onesys_Y(:,3)=onesys_Y(:,3)+1000;  %re-order market ID to avoid replication of MKT ID with the new market ID
       baseY_sys=[baseY_sys;onesys_Y];        
    end
    % construct X for affiliated hospitals' brothers outside the selected markets
    onesys_X=tempX(tempX(:,end)==sysid(i),:);
    incmkt=onesys_X(:,3);  %columnn vector
    exclu=max((~(incmkt(:,ones(numel(thisboot),1))-thisboot(ones(size(onesys_X,1),1),:))),[],2);  %only include bros OUTSIDE the selected markets
    onesys_X(exclu==1,:)=[];
    if isempty(onesys_X)==0
       onesys_X(:,3)=onesys_X(:,3)+1000;  %re-order market ID to avoid replication of MKT ID with the new market ID        
       baseX_sys=[baseX_sys;onesys_X];        
    end
    % construct Z for affiliated hospitals' brothers outside the selected markets
    onesys_Z=tempZ(tempZ(:,end)==sysid(i),:);
    incmkt=onesys_Z(:,4);  %columnn vector
    exclu=max((~(incmkt(:,ones(numel(thisboot),1))-thisboot(ones(size(onesys_Z,1),1),:))),[],2);  %only include bros OUTSIDE the selected markets
    onesys_Z(exclu==1,:)=[];
    if isempty(onesys_Z)==0
       onesys_Z(:,4)=onesys_Z(:,4)+1000;  %re-order market ID to avoid replication of MKT ID with the new market ID                
       baseZ_sys=[baseZ_sys;onesys_Z];        
    end
end

% attach outside bros
baseY=[baseY;baseY_sys];
baseX=[baseX;baseX_sys];
baseZ=[baseZ;baseZ_sys];
    
% create new ahaid (first sort baseY, baseX and baseZ to make sure they are
% in the SAME order)
baseY=sortrows(baseY,[3,1,2]);
baseX=sortrows(baseX,[3,1,2]);
baseZ=sortrows(baseZ,[4,1,2,3]);
[~,~,newid]=unique(baseY(:,[3,1]),'rows');  % Creating IDs by groups in the order of [hsa, ahaid]
corr_newid_ahaid_newhsa=unique([newid,baseY(:,[1,3])],'rows');  % Crosswork between newid, ahaid, and reorder_mkt
id_baseZ=kron(newid,ones(ns,1));

% plug in new ahaid
baseY=[newid,baseY(:,2:end)];
baseX=[newid,baseX(:,2:end)];
baseZ=[id_baseZ,baseZ(:,2:end)];


% ### UPDATED on 12/03/2020: CANNOT use the vmshbd1-4 and vlag1-4 because {bed25,
% bed50, bed75} vary by bootstrap samples (thisboot)
bed25=prctile(baseX(baseX(:,end)==0,4),25);
bed50=median(baseX(baseX(:,end)==0,4));
bed75=prctile(baseX(baseX(:,end)==0,4),75);

all_bdtype=[(baseX(:,4)<=bed25&baseX(:,4)>0), (baseX(:,4)<=bed50&baseX(:,4)>bed25), (baseX(:,4)<=bed75&baseX(:,4)>bed50), (baseX(:,4)>bed75)];
baseZ(:,5:8)=repmat(baseZ(:,5)+baseZ(:,6)+baseZ(:,7)+baseZ(:,8),1,4).*kron(all_bdtype,ones(ns,1));
baseZ(:,9:12)=repmat(baseZ(:,9)+baseZ(:,10)+baseZ(:,11)+baseZ(:,12),1,4).*kron(all_bdtype,ones(ns,1));


% ### UPDATED on 12/03/2020: Re-calculate vsys because affiliated hospitals
% in one system could vary by bootstrap samples (thisboot)
sys_in_thisboot=unique(baseY(:,end));
sys_in_thisboot(sys_in_thisboot==0)=[];
update_vsysh=zeros(size(baseY,1),ns);
for year=2006:2009
    for i=1:numel(sys_in_thisboot)
        tmpsys=baseY(baseY(:,end)==sys_in_thisboot(i)&baseY(:,2)==year,5);  
        nummmb=numel(tmpsys);
        tmpsys(tmpsys==1)=[];                
        if isempty(tmpsys)==0
           occur=histc(tmpsys,choice);
           occur=reshape(occur,1,[]); 
           onevsysh=occur/sum(occur);
           update_vsysh(baseY(:,end)==sys_in_thisboot(i)&baseY(:,2)==year,:)=onevsysh((1:size(onevsysh,1))' * ones(1,nummmb), (1:size(onevsysh,2))' );
        end    
    end
end
baseZ(:,14)=reshape(update_vsysh',[],1);


%%%$$$%%% create inputs for newmain file (structural estimation)
% ### UPDATED on 12/06/2020: generage "regvar6_sing.raw"-"regvar9_sing.raw"
baseY=sortrows(baseY,[1,2]);
baseX=sortrows(baseX,[1,2]);

regvar=[baseX(baseX(:,end)==0,[1,2,4,5]),baseY(baseY(:,end)==0,[4,3])];   %matrix: newahaid year bdtot prechs adj_vnd hsa
for i=6:9
    textFileName_regvar=sprintf('regvar%d_sing.raw',i);
    output_regvar=regvar(regvar(:,2)==2000+i,:);
    output_regvar=sortrows(output_regvar,1);
    dlmwrite(textFileName_regvar,output_regvar,'delimiter', ',', 'precision', 8)
end

%the matrix: ahaid year hhi mktcat totpop65
ind_mktcat=baseX(:,13:15);
mktcat=zeros(size(baseX,1),1);
for i=1:size(baseX,1)
    if sum(ind_mktcat(i,:))>0
        mktcat(i)=find(ind_mktcat(i,:)>0);        
    end    
end
textFileName_mktcat=sprintf('mktvar_pophhi.raw');
dlmwrite(textFileName_mktcat,[baseX(baseX(:,end)==0,[1,2,9]),mktcat(baseX(:,end)==0,:),baseX(baseX(:,end)==0,10)],'delimiter', ',', 'precision', 8)
%%%$$$%%% END



% create inputs for others; ### UPDATED on 06/08/2020: only keep bdtot for
% non-adopting hospitals to avoid non-concavity
Y_SANA=baseY(baseY(:,end)==0&baseY(:,end-1)==1,:);
Y_SANA=Y_SANA(:,[1,2,4]);  %matrix: [newahaid year adj_vnd]
X_SANA=baseX(baseX(:,end)==0&baseX(:,5)==1,:);
%X_SANA=X_SANA(:,[1,2,4,7,8,9,10]);   %matrix: [newahaid year bdtot profit teaching hhi totpop65]; ### changed from nprofit to profit to avoid non-concavity
X_SANA=X_SANA(:,[1,2,4,9,10,16]);   %matrix: [newahaid year bdtot hhi totpop65 yrsince]; ### only keep "bdtot" to avoid non-concavity
Z_SANA=baseZ(baseZ(:,end)==0&baseZ(:,end-1)==1,:);
Z_SANA=Z_SANA(:,[1,2,3,5:8,15:17,18]);  %matrix: [newahaid year vid vmshbd1 vmshbd2 vmshbd3 vmshbd4 below50cat below75cat below100cat yrsinceXDOM]

Y_SAAD=baseY(baseY(:,end)==0&baseY(:,end-1)>1,:);
Y_SAAD=Y_SAAD(:,[1,2,4]); %matrix: [newahaid year adj_vnd]
X_SAAD=baseX(baseX(:,end)==0&baseX(:,5)>1,:); 
X_SAAD=X_SAAD(:,[1,2,4,7,8,9,10,16]);   %matrix: [newahaid year bdtot profit teaching hhi totpop65 yrsince]
Z_SAAD=baseZ(baseZ(:,end)==0&baseZ(:,end-1)>1,:);
Z_SAAD=Z_SAAD(:,[1,2,3,5:8,9:12,18]);  %matrix: [newahaid year vid vmshbd1 vmshbd2 vmshbd3 vmshbd4 vlagbd1 vlagbd2 vlagbd3 vlagbd4 yrsinceXDOM]

Y_AFNA=baseY(baseY(:,end)>0&baseY(:,end-1)==1,:);
Y_AFNA=Y_AFNA(:,[1,2,4]); %matrix: [newahaid year adj_vnd]
X_AFNA=baseX(baseX(:,end)>0&baseX(:,5)==1,:);
%X_AFNA=X_AFNA(:,[1,2,4,6,11,9,10]);   %matrix: [newahaid year bdtot nprofit perc_mcr hhi totpop65]
X_AFNA=X_AFNA(:,[1,2,4,9,10,16]);   %only keep "bdtot"; matrix: [newahaid year bdtot hhi totpop65 yrsince]
Z_AFNA=baseZ(baseZ(:,end)>0&baseZ(:,end-1)==1,:);
Z_AFNA=Z_AFNA(:,[1,2,3,14,15:17]);  %matrix: [newahaid year vid newvsysh below50cat below75cat below100cat]

Y_AFAD=baseY(baseY(:,end)>0&baseY(:,end-1)>1,:);
Y_AFAD=Y_AFAD(:,[1,2,4]); %matrix: [newahaid year adj_vnd]
X_AFAD=baseX(baseX(:,end)>0&baseX(:,5)>1,:);
X_AFAD=X_AFAD(:,[1,2,4,11,12,9,10,16]);   %matrix: [newahaid year bdtot perc_mcr perc_mcd hhi totpop65 yrsince]
Z_AFAD=baseZ(baseZ(:,end)>0&baseZ(:,end-1)>1,:);
Z_AFAD=Z_AFAD(:,[1,2,3,14,13]);  %matrix: [newahaid year vid newvsysh newvlag]




%&&& Obtain estimates from external functions: clogit and pclogit
%!!! for affiliated hospitals: remember to remove those out of the
%considered boot!!!

%&&&&&&&&&& SANA
b_SANA=Clogit_SANA_decsunk(Y_SANA,X_SANA,Z_SANA);
P_SANA=Pclogit_SANA_decsunk(Y_SANA,X_SANA,Z_SANA,b_SANA);

%test coefficients
%{
textFileNameY_sana=sprintf('testYsana.csv');
textFileNameX_sana=sprintf('testXsana.csv');
textFileNameZ_sana=sprintf('testZsana.csv');

dlmwrite(textFileNameY_sana,Y_SANA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameX_sana,X_SANA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameZ_sana,Z_SANA,'delimiter', ',', 'precision', 8)
%}

%produce coefficients for simulation
nb = size(X_SANA,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_SANA,2)-3;
nc=numel(unique(Y_SANA(:,end)));

alpha=[b_SANA(nb*(nc-1)+1:nb*(nc-1)+4)' 0 0 0 0 b_SANA(end)];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q
mktcat_alpha=[b_SANA(end-3) b_SANA(end-2) b_SANA(end-1)];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_alpha=[zeros(numel(mktcat_alpha),1),repmat(mktcat_alpha',1,ns-1)];
xb=reshape(b_SANA(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

gamma2=xb(1,:);
gamma3=xb(2,:);
gamma4=xb(3,:);
gamma5=xb(4,:);
gamma6=xb(5,:);
gamma7=xb(6,:);
gamma8=xb(7,:);
gamma9=xb(8,:);
gamma10=xb(9,:);
gamma11=xb(10,:);
gamma12=xb(11,:);
gamma13=xb(12,:);

Gamma=[zeros(size(gamma2))' gamma2' gamma3' gamma4' gamma5' gamma6' gamma7' gamma8' gamma9' gamma10' gamma11' gamma12' gamma13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)
clear xb nb nb2 nc

hosp_SANA=baseX(baseX(:,end)==0&baseX(:,5)==1,:);
hosp_SANA=[hosp_SANA(:,1:5),P_SANA,hosp_SANA(:,[6:end-2,end])];  %### UPDATED on 12/03/2020: Skip yrsince



%&&&&&&&&&& SAAD
b_SAAD=Clogit_SAAD_decsunk(Y_SAAD,X_SAAD,Z_SAAD);
P_SAAD=Pclogit_SAAD_decsunk(Y_SAAD,X_SAAD,Z_SAAD,b_SAAD);


%test coefficients
%{
textFileNameY_saad=sprintf('testYsaad.csv');
textFileNameX_saad=sprintf('testXsaad.csv');
textFileNameZ_saad=sprintf('testZsaad.csv');

dlmwrite(textFileNameY_saad,Y_SAAD,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameX_saad,X_SAAD,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameZ_saad,Z_SAAD,'delimiter', ',', 'precision', 8)
%}


%produce coefficients for simulation
nb = size(X_SAAD,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_SAAD,2)-3;
nc=numel(unique(Y_SAAD(:,end)));

tau=b_SAAD(nb*(nc-1)+1:nb*(nc-1)+9)' ; %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q
mktcat_tau=[0	0	0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

xb=reshape(b_SAAD(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

zeta2=[0 0 0 0 0 0 0];
zeta3=xb(1,:);
zeta4=xb(2,:);
zeta5=xb(3,:);
zeta6=xb(4,:);
zeta7=xb(5,:);
zeta8=xb(6,:);
zeta9=xb(7,:);
zeta10=xb(8,:);
zeta11=xb(9,:);
zeta12=xb(10,:);
zeta13=xb(11,:);

Zeta=[zeros(size(zeta2))' zeta2' zeta3' zeta4' zeta5' zeta6' zeta7' zeta8' zeta9' zeta10' zeta11' zeta12' zeta13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)
clear xb nb nb2 nc

hosp_SAAD=baseX(baseX(:,end)==0&baseX(:,5)>1,:);
hosp_SAAD=[hosp_SAAD(:,1:5),[zeros(size(P_SAAD,1),1),P_SAAD],hosp_SAAD(:,[6:end-2,end])]; %### UPDATED on 12/03/2020: Skip yrsince
 


%&&&&&&&&&& AFNA 
b_AFNA=Clogit_AFNA_decsunk(Y_AFNA,X_AFNA,Z_AFNA);
P_AFNA=Pclogit_AFNA_decsunk(Y_AFNA,X_AFNA,Z_AFNA,b_AFNA);


%test coefficients
%{
textFileNameY_afna=sprintf('testYafna.csv');
textFileNameX_afna=sprintf('testXafna.csv');
textFileNameZ_afna=sprintf('testZafna.csv');

dlmwrite(textFileNameY_afna,Y_AFNA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameX_afna,X_AFNA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameZ_afna,Z_AFNA,'delimiter', ',', 'precision', 8)
%}

%produce coefficients for simulation
nb = size(X_AFNA,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_AFNA,2)-3;
nc=numel(unique(Y_AFNA(:,end)));

sigma=[0 0 0 b_AFNA(nb*(nc-1)+1)];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q
mktcat_sigma=[b_AFNA(end-2) b_AFNA(end-1) b_AFNA(end)];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_sigma=[zeros(numel(mktcat_sigma),1),repmat(mktcat_sigma',1,ns-1)];

xb=reshape(b_AFNA(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

eta2=xb(1,:);
eta3=xb(2,:);
eta4=xb(3,:);
eta5=xb(4,:);
eta6=xb(5,:);
eta7=xb(6,:);
eta8=xb(7,:);
eta9=xb(8,:);
eta10=xb(9,:);
eta11=xb(10,:);
eta12=xb(11,:);
eta13=xb(12,:);

Eta=[zeros(size(eta2))' eta2' eta3' eta4' eta5' eta6' eta7' eta8' eta9' eta10' eta11' eta12' eta13'];

hosp_AFNA=baseX(baseX(:,end)>0&baseX(:,5)==1,:);
hosp_AFNA=[hosp_AFNA(:,1:5),P_AFNA,hosp_AFNA(:,[6:end-2,end])]; %### UPDATED on 12/03/2020: Skip yrsince



%&&&&&&&&&& AFAD 
b_AFAD=Clogit_AFAD_decsunk(Y_AFAD,X_AFAD,Z_AFAD);
P_AFAD=Pclogit_AFAD_decsunk(Y_AFAD,X_AFAD,Z_AFAD,b_AFAD);


%test coefficients
%{
textFileNameY_afad=sprintf('testYafad.csv');
textFileNameX_afad=sprintf('testXafad.csv');
textFileNameZ_afad=sprintf('testZafad.csv');

dlmwrite(textFileNameY_afad,Y_AFAD,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameX_afad,X_AFAD,'delimiter', ',', 'precision', 8)
dlmwrite(textFileNameZ_afad,Z_AFAD,'delimiter', ',', 'precision', 8)
%}

%produce coefficients for simulation
nb = size(X_AFAD,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_AFAD,2)-3;
nc=numel(unique(Y_AFAD(:,end)));

rho=[0 0 b_AFAD(end) b_AFAD(end-1)]; %the coefficient: vmkt vmsh vlag vsys
mktcat_rho=[0    0   0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

xb=reshape(b_AFAD(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

lambda2=[0 0 0 0 0 0 0]; %the order is as follows: preadp(indicate whether adopting or not) bdtot perc_mcr perc_mcd hhi totpop65 constant
lambda3=xb(1,:);
lambda4=xb(2,:);
lambda5=xb(3,:);
lambda6=xb(4,:);
lambda7=xb(5,:);
lambda8=xb(6,:);
lambda9=xb(7,:);
lambda10=xb(8,:);
lambda11=xb(9,:);
lambda12=xb(10,:);
lambda13=xb(11,:);

Lambda=[zeros(size(lambda2))' lambda2' lambda3' lambda4' lambda5' lambda6' lambda7' lambda8' lambda9' lambda10' lambda11' lambda12' lambda13'];

hosp_AFAD=baseX(baseX(:,end)>0&baseX(:,5)>1,:);
hosp_AFAD=[hosp_AFAD(:,1:5),[zeros(size(P_AFAD,1),1),P_AFAD],hosp_AFAD(:,[6:end-2,end])]; %### UPDATED on 12/03/2020: Skip yrsince




%For debugging: EXPORT parameters
%{
textFileName_para_SANA=sprintf('para_SANA.csv');
textFileName_para_SAAD=sprintf('para_SAAD.csv');
textFileName_para_AFNA=sprintf('para_AFNA.csv');
textFileName_para_AFAD=sprintf('para_AFAD.csv');

textFileName_hosp_SANA=sprintf('hosp_SANA.csv');
textFileName_hosp_SAAD=sprintf('hosp_SAAD.csv');
textFileName_hosp_AFNA=sprintf('hosp_AFNA.csv');
textFileName_hosp_AFAD=sprintf('hosp_AFAD.csv');

dlmwrite(textFileName_para_SANA,b_SANA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName_para_SAAD,b_SAAD,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName_para_AFNA,b_AFNA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName_para_AFAD,b_AFAD,'delimiter', ',', 'precision', 8)

dlmwrite(textFileName_hosp_SANA,hosp_SANA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName_hosp_SAAD,hosp_SAAD,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName_hosp_AFNA,hosp_AFNA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName_hosp_AFAD,hosp_AFAD,'delimiter', ',', 'precision', 8)



%For debugging: IMPORT the saved parameters and hospital info

textFileName_para_SANA=sprintf('para_SANA.csv');
textFileName_para_SAAD=sprintf('para_SAAD.csv');
textFileName_para_AFNA=sprintf('para_AFNA.csv');
textFileName_para_AFAD=sprintf('para_AFAD.csv');

textFileName_hosp_SANA=sprintf('hosp_SANA.csv');
textFileName_hosp_SAAD=sprintf('hosp_SAAD.csv');
textFileName_hosp_AFNA=sprintf('hosp_AFNA.csv');
textFileName_hosp_AFAD=sprintf('hosp_AFAD.csv');

b_SANA=dlmread(textFileName_para_SANA);
nb = size(X_SANA,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_SANA,2)-3;
nc=numel(unique(Y_SANA(:,end)));

alpha=[b_SANA(nb*(nc-1)+1:nb*(nc-1)+4)' 0 0 0 0 b_SANA(end)];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q
mktcat_alpha=[b_SANA(end-3) b_SANA(end-2) b_SANA(end-1)];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_alpha=[zeros(numel(mktcat_alpha),1),repmat(mktcat_alpha',1,ns-1)];
xb=reshape(b_SANA(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

gamma2=xb(1,:);
gamma3=xb(2,:);
gamma4=xb(3,:);
gamma5=xb(4,:);
gamma6=xb(5,:);
gamma7=xb(6,:);
gamma8=xb(7,:);
gamma9=xb(8,:);
gamma10=xb(9,:);
gamma11=xb(10,:);
gamma12=xb(11,:);
gamma13=xb(12,:);

Gamma=[zeros(size(gamma2))' gamma2' gamma3' gamma4' gamma5' gamma6' gamma7' gamma8' gamma9' gamma10' gamma11' gamma12' gamma13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)
clear xb nb nb2 nc



b_SAAD=dlmread(textFileName_para_SAAD);
nb = size(X_SAAD,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_SAAD,2)-3;
nc=numel(unique(Y_SAAD(:,end)));


tau=b_SAAD(nb*(nc-1)+1:nb*(nc-1)+9)' ; %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q
mktcat_tau=[0	0	0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

xb=reshape(b_SAAD(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

zeta2=[0 0 0 0 0 0 0];
zeta3=xb(1,:);
zeta4=xb(2,:);
zeta5=xb(3,:);
zeta6=xb(4,:);
zeta7=xb(5,:);
zeta8=xb(6,:);
zeta9=xb(7,:);
zeta10=xb(8,:);
zeta11=xb(9,:);
zeta12=xb(10,:);
zeta13=xb(11,:);

Zeta=[zeros(size(zeta2))' zeta2' zeta3' zeta4' zeta5' zeta6' zeta7' zeta8' zeta9' zeta10' zeta11' zeta12' zeta13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)
clear xb nb nb2 nc


b_AFNA=dlmread(textFileName_para_AFNA);
nb = size(X_AFNA,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_AFNA,2)-3;
nc=numel(unique(Y_AFNA(:,end)));

sigma=[0 0 0 b_AFNA(nb*(nc-1)+1)];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q
mktcat_sigma=[b_AFNA(end-2) b_AFNA(end-1) b_AFNA(end)];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_sigma=[zeros(numel(mktcat_sigma),1),repmat(mktcat_sigma',1,ns-1)];

xb=reshape(b_AFNA(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

eta2=xb(1,:);
eta3=xb(2,:);
eta4=xb(3,:);
eta5=xb(4,:);
eta6=xb(5,:);
eta7=xb(6,:);
eta8=xb(7,:);
eta9=xb(8,:);
eta10=xb(9,:);
eta11=xb(10,:);
eta12=xb(11,:);
eta13=xb(12,:);

Eta=[zeros(size(eta2))' eta2' eta3' eta4' eta5' eta6' eta7' eta8' eta9' eta10' eta11' eta12' eta13'];


b_AFAD=dlmread(textFileName_para_AFAD);
nb = size(X_AFAD,2)-2+1;  %add the constant, in addition to the covariates
nb2= size(Z_AFAD,2)-3;
nc=numel(unique(Y_AFAD(:,end)));

rho=[0 0 b_AFAD(end) b_AFAD(end-1)]; %the coefficient: vmkt vmsh vlag vsys
mktcat_rho=[0    0   0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

xb=reshape(b_AFAD(1:nb*(nc-1)),nb,[])';
xb=[xb(:,2:end),xb(:,1)];

lambda2=[0 0 0 0 0 0 0]; %the order is as follows: preadp(indicate whether adopting or not) bdtot perc_mcr perc_mcd hhi totpop65 constant
lambda3=xb(1,:);
lambda4=xb(2,:);
lambda5=xb(3,:);
lambda6=xb(4,:);
lambda7=xb(5,:);
lambda8=xb(6,:);
lambda9=xb(7,:);
lambda10=xb(8,:);
lambda11=xb(9,:);
lambda12=xb(10,:);
lambda13=xb(11,:);
Lambda=[zeros(size(lambda2))' lambda2' lambda3' lambda4' lambda5' lambda6' lambda7' lambda8' lambda9' lambda10' lambda11' lambda12' lambda13'];

hosp_SANA=dlmread(textFileName_hosp_SANA);
hosp_SAAD=dlmread(textFileName_hosp_SAAD);
hosp_AFNA=dlmread(textFileName_hosp_AFNA);
hosp_AFAD=dlmread(textFileName_hosp_AFAD);
%}



%&&&&&&&&&& Start simulation
%bed25=25;
%bed50=80;
%bed75=193;
%bed25=prctile(baseX(baseX(:,end)==0,4),25);
%bed50=median(baseX(baseX(:,end)==0,4));
%bed75=prctile(baseX(baseX(:,end)==0,4),75);
%ns=1;       %will take 516 secs if ns=500
T=100;
%s=1;
textFileName3=sprintf('fs%d.csv',s);  %the output file
textFileName4=sprintf('ddom%d.csv',s);   %the output file
textFileName5=sprintf('eerr%d.csv',s);   %the output file
textFileName6=sprintf('mks%d.csv',s);   %the output file

rng(s);
FSTA=[];
DDOM=[];
EERR=[]; %### UPDATED on 03/16/2017: the output including the probability as the expected error
MKTS=[]; %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation

for k=2006:2009
%k=7;
P=[hosp_SANA;hosp_SAAD];
P(P(:,2)~=k,:)=[];
P(:,2)=[];
PP=[hosp_AFNA;hosp_AFAD];
PP(PP(:,2)~=k,:)=[];
PP(:,2)=[];
allsys=PP(:,end);
%The following bracket section may be DELETED
%{
textFileName1a=sprintf('prdprob_a%d.csv',k); %the input file
textFileName1na=sprintf('prdprob_na%d.csv',k); %the input file
textFileName2a=sprintf('syschs_a%d.csv',k);  %the input file
textFileName2na=sprintf('syschs_na%d.csv',k);  %the input file
J1=dlmread(textFileName1na);  % ### UPDATED on 04/21/2020---the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName1a); 
%reshape the raw data from STATA, SA-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
H1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, SA-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
H2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
P=[H1;H2];
%pick=find(P(:,2)==17 | P(:,2)==10 |P(:,2)==44 |P(:,2)==46 |P(:,2)==62 |P(:,2)==90);
%P=P(pick,:);
J1=dlmread(textFileName2na); %the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName2a); 
%reshape the raw data from STATA, AF-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
HH1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, AF-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
HH2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
PP=[HH1;HH2];
allsys=PP(:,end);
%}
%Specify each page of FSTA and DDOM: fs and G
fs=zeros(size(P,1)*size(choice,2),T+5); 
G=zeros(size(P,1)*size(choice,2),T+4);
MS=zeros(size(P,1)*size(choice,2),T+4); %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation
eer=zeros(size(P,1)*ns,T+3);  %### UPDATED on 03/16/2017: include the simulated expected errors: ahaid, year, option given, prob@t=1, prob@t=2, ..., prob@t=100

for i=1:size(P,1) 
    tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
    exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
    exomkt=sortrows(exomkt,1);
    altemp=[exomkt;tempHH];  %group all the rival hospitals (standalone and affiliated) together (to do the random draw)
    tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
    alltemp=[tempmkt;tempHH]; %group all the hospitals (including the one considered) (standalone and affiliated) together (to do the random draw)
    %The following is to compute fs0 and g0. Then add [year, bdtot, fs0]
    %into fs and add [year,g0] into G.
    fs0=P(P(:,1)==P(i,1),4);
    mkts0=0;
    prechs=altemp(:,4);
    prechs(prechs==1)=[];
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
       mktsh=occur/sum(occur);
       mkts0=mktsh(fs0);
    end    
    prechs(prechs==2)=[];
    prechs(prechs==13)=[];
    if isempty(prechs)==1
       g0=0;
    end
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
     %  if occur(2)>0
     %     occur(2)=1;
     %  end
       occur=reshape(occur,1,[]);
       mktdom=choice(occur==max(occur));
       g0=max(fs0-mktdom==0);
    end
  %%%%%%%%%%%%%%%%%%%%%%%%%
  %consider the case where the hospital has already installed EMR  
   if P(i,4)>1   
        want=zeros(size(choice,2)-1,T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2)-1,T+1);  
        mkts=zeros(size(choice,2)-1,T+1);
        tempeer=zeros(ns-1,T+1);  %### UPDATED on 03/16/2017: a section of eer: the probability at each period
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market. ### 04/21/2020: should be
        %for stand-alone hospitals whose competitors are ALL stand-alone
        %hospitals, too.
        if isempty(sys)==1
          for l=2:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;
              if isempty(altemp)==1
                  newchs=l;
              end
              if isempty(altemp)==0
                 w=altemp(:,5:17); 
                 cc=cumsum(w,2);
                 cc(:,end)=1;
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1  %for monopoly market
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l;
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END                      
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                                        
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0;
                       %mkts(l-1,t+1)=0;
                     end
                     if isempty(prechs)==0
                        occur=histc(prechs,choice)';
                        occur=reshape(occur,1,[]);           
                        vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                        %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                        %vmshf=vmshff;
                        prechs(prechs==2)=[];
                        prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0; %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4);
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);  %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".                   
                     Tchoice=choice';                     
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);
                    
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';    
                     % ### UPDATED on 12/04/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t-(k-t0))*(t+k-t0<=tbar)+zeros(size(vmktf,1),1)*(t+k-t0>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;  %only for DEBUG
                     % ### UPDATED on 06/08/2020: for all non-adopting
                     % hospitals, (SANA and AFNA), only keep bdtot and
                     % market variables (HHI and totpop65)
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,3)*Gamma(1,:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([end-3,end-2],:)... 
                       +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(end-1,:);                        
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:)...
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);                      
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     %v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                      % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Gamma([1,2,3],:)...
                      % +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:);                  
                     %v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                       %  +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                        % +alltemp(:,[21,22])*Zeta([4,5],:);                         
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                                        
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);  %###!!!If the first one changed back to using max, remember to use "altemp"!
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];                   
                     %### UPDATED on 1/4/2017: re-define "selfnewchs" by
                     %figuring out what brings in largest payoff given the
                     %rivals' new choice, from "newchs(2:end)".            
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),k*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),k*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),k*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),k*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market. ### 04/21/2020: consider hospitals whose competitors
        %include affiliated hospitals
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:); %find out the outside ones
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=2:13
              %w=altemp(:,5:17);
              %chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1;
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l; 
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END                      
                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     
                    alltemp(:,4)=newchs(1:size(alltemp,1),:);
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                     
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0; %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0;  %### UPDATED on 07/16/2018                     
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                  
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                    vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                    vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                    

                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);                           
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);                    
                     %bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);                     
                     %bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)>1);                                      
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                 
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                     % +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                     % +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                        
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                    
                     % ### UPDATED on 12/04/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t-(k-t0))*(t+k-t0<=tbar)+zeros(size(vmktf,1),1)*(t+k-t0>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;  %only for DEBUG                     
                     % ### UPDATED on 06/08/2020: for all non-adopting
                     % hospitals, (SANA and AFNA), only keep bdtot and
                     % market variables (HHI and totpop65)
                      v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                      +cons1(ones(size(vmktf,1),1),:)+univ(:,3)*Gamma(1,:)...
                      +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([end-3,end-2],:)...
                      +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(end-1,:);                                                        
                      v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([end-3,end-2],:)... 
                      +yrsince*Eta(end-1,:);                 
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary 
                     %v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                      % +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Gamma([1,2,3],:)...
                      % +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:)...
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);                                         
                    % v3=vsysh*sigma(4)...
                    %  +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                    %  +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:);                  
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:)...
                      +yrsince*Lambda(6,:);
                  
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                     
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);  %!!! change it to "univ" or "exouniv" depends on random or max.
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),k*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),k*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),k*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),k*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
   end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %consider the case where the hospital has never installed EMR
   if P(i,4)==1  
        %tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
        %exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
        %exomkt=sortrows(exomkt,1);
        %altemp=[exomkt;tempHH];  %group all the hospitals (standalone and affiliated) together to do the random draw)
        %tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pc1 pc2...pc13
        %alltemp=[tempmkt;tempHH];
        want=zeros(size(choice,2),T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2),T+1); 
        mkts=zeros(size(choice,2),T+1);
        tempeer=zeros(ns,T+1); 
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market.
        if isempty(sys)==1
          for l=1:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;             
              if isempty(altemp)==1  %consider the case of monopoly
                  newchs=l;
              end
              if isempty(altemp)==0
                 w=altemp(:,5:17); 
                 cc=cumsum(w,2);
                 cc(:,end)=1;                 
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];
                                          
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;   %### UPDATED on 07/16/2018
                     end
                     
                     
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);     %### UPDATED on 07/16/2018                 
                     Tchoice=choice';      
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                   
                   curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);


                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);    
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);                  
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);                  
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';   
                     % ### UPDATED on 12/04/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t-(k-t0))*(t+k-t0<=tbar)+zeros(size(vmktf,1),1)*(t+k-t0>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;  %only for DEBUG                     
                     % ### UPDATED on 06/08/2020: for all non-adopting
                     % hospitals, (SANA and AFNA), only keep bdtot and
                     % market variables (HHI and totpop65)
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,3)*Gamma(1,:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([end-3,end-2],:)... 
                       +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(end-1,:);                        
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:)...
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);                                                           
                      % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     %v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                      % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Gamma([1,2,3],:)...
                      % +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:);                  
                   %  v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                    %     +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                    %     +alltemp(:,[21,22])*Zeta([4,5],:);                 
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;              
                     %selfnewchs=find(w(1,:)==max(w(1,:)));
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l,t+1)=newchs(1);  
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),k*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),k*ones(size(g,1),1),g0*ones(size(g,1),1)],g]; 
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),k*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];  %### UPDATED on 07/16/2018
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),k*ones(size(tempeer,1),1),tempeer];
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market.
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that have belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:);
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=1:13
              %w=altemp(:,5:17);
              %chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1; 
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END

                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     alltemp(:,4)=newchs(1:size(alltemp,1),:);
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];                               
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;   %### UPDATED on 07/16/2018                   
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                 
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                     vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                     vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                   
                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);  
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);    
                                     
                    % bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);     
                    % bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)==1);                                        
                    % v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                    %  +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                   
                    % v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                    %  +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                    % v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                    %  +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                    % v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                    %  +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                    
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';
                     % ### UPDATED on 12/04/2020: generate yrsince
                     yrsince=ones(size(vmktf,1),1)*(tbar-t-(k-t0))*(t+k-t0<=tbar)+zeros(size(vmktf,1),1)*(t+k-t0>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;  %only for DEBUG 
                     % ### UPDATED on 06/08/2020: for all non-adopting
                     % hospitals, (SANA and AFNA), only keep bdtot and
                     % market variables (HHI and totpop65)
                      v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                      +cons1(ones(size(vmktf,1),1),:)+univ(:,3)*Gamma(1,:)...
                      +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([end-3,end-2],:)...
                      +vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(end-1,:);                                                        
                      v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([end-3,end-2],:)... 
                      +yrsince*Eta(end-1,:);                 
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary 
                     %v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                      % +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Gamma([1,2,3],:)...
                      % +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:)...
                         +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:);                                         
                    % v3=vsysh*sigma(4)...
                    %  +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                    %  +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:);                  
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:)...
                      +yrsince*Lambda(6,:);
                    
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;
                     %{
                     selfnewchs=find(w(1,:)==max(w(1,:)));
                     w(1,:)=[];
                     cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l,t+1)=newchs(1); 
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),k*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),k*ones(size(g,1),1),g0*ones(size(g,1),1)],g];   
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),k*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),k*ones(size(tempeer,1),1),tempeer];

        end
    end
end
G(:,[3,4])=G(:,[4,3]);  %exchange g0 with action_given
MS(:,[3,4])=MS(:,[4,3]); % ### UPDATED on 08/22/2018: exchange mkts0 with action_given
FSTA=[FSTA;fs];
DDOM=[DDOM;G];
MKTS=[MKTS;MS];
EERR=[EERR;eer];
end
dlmwrite(textFileName3,FSTA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName4,DDOM,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName5,EERR,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName6,MKTS,'delimiter', ',', 'precision', 8)


toc


